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Summary 

This article focuses on construction, directly in physical space, of three-dimensional covariance 
functions parametrized by a tunable length field, and on an application of this theory to reproduce the 
Quasi-Biennial Oscillation (QBO) in the Goddard Earth Observing System, Version 4 (GEOS-4) data 
assimilation system. These covariance models are referred to as multi-level or nonseparable, to associate 
them with the application where a multi-level covariance with a large troposphere to stratosphere length 
field gradient is used to reproduce the QBO from sparse radiosonde observations in the tropical lower 
stratosphere. The multi-level covariance functions extend well-known single-level covariance functions 
depending only on a length scale. Generalizations of the first- and third-order autoregressive covariances 
in three dimensions are given, providing multi-level covariances with zero and three derivatives at zero 
separation, respectively. Multi-level piecewise rational covariances with two continuous derivatives at zero 
separation are also provided. Multi-level powerlaw covariances are constructed with continuous deriva- 
tives of all orders. Additional multi-level covariance functions are constructed using the Schur product 
of single- and multi-level covariance functions. A multi-level powerlaw covariance used to reproduce the 
QBO in GEOS-4 is described along with details of the assimilation experiments. The new covariance 
model is shown to represent the vertical wind shear associated with the QBO much more effectively than 
in the baseline GEOS-4 system. 
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1 Introduction 

This article focuses on constructive development and application of covariance func- 
tions with length scales determined by certain prescribed functions. These reduce to 
covariance functions, depending only on a length scale, that have been widely applied 
in diverse applications where covariance models are required. Covariance modeling ap- 
plications in the earth sciences are emphasized, where correlation length scales vary in 
altitude and/or on horizontal surfaces. In accordance with standard terminology, these 
covariance functions will be said to depend on a variable length field (cf. Yaglom 1987, 
p. 19). 

Parametrized forecast- and observation-error covariance function models depending 

on a length field have been an important component of atmospheric data assimilation 
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systems for many years. Efforts to construct covariance functions with vertically variable 
length fields evolved in response to empirical studies of short-range forecast error statistics 
by Hollingsworth and Lonnberg (1986) and Lonnberg and Hollingsworth (1986), where 
it was shown that horizontal forecast error correlation length scales increase significantly 
with height. The terminology nonseparable is used throughout the literature to describe 
covariance models with length fields that vary with height (Phillips 1986; Bartello and 
Mitchell 1992; Courtier et al. 1998; Rabier et al. 1998; Andersson et al. 1998; Daley and 
Barker 2000), since three-dimensional covariance functions with independent horizon- 
tal and vertical components are called separable (cf. Daley 1991, Sec. 4.3). Short-range 
forecast errors are known to be smaller over data-dense areas like the northern hemi- 
sphere continents than over data-poor areas like the oceans (Bouttier 1994; Courtier et 
al 1998, Sec. 3(vi); Andersson et al 1998). The covariance functions described below 
could be applied, for example, with short correlation lengths over data-dense areas and 
long correlation lengths over the oceans. 

Forecast and observation error covariance models in 3D-Var systems are fixed (flow 
independent) over the entire assimilation, and subject to one or more simplifying as- 
sumptions such as homogeneity, isotropy, compact support, and separability (Parrish 
and Derber 1992; Cohn et al. 1998; Courtier et al. 1998; Gaspari and Cohn 1999; Daley 
and Barker 2000). These limitations are also present in 4D-Var systems, since the initial 
covariance at each analysis cycle is specified as in 3D-Var (Rabier et al. 2000; Mahfouf et 
al 2000; Klinker et aL 2000). The new covariance functions constructed below directly 
extend the single-level covariance functions within the 3D-Var like Physical-space Sta- 
tistical Analysis System (PSAS; Cohn et al. 1998). Since the new covariance functions 
extend models that are widely employed, development efforts needed to incorporate them 
into other existing 3D-Var and 4D-Var systems should be modest. 
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Results are presented in two sections. Section 2 develops mathematical results needed 
for the construction of global, three-dimensional covariance functions directly in physical 
space, that depend on a length field. The style of presentation mirrors Caspar! and Cohn 
(1999, henceforth referred to as GC), with extensive use of results from that article. A 
first reading of GC is highly recommended as background motivating the results in Sec. 2. 
Although the theory and constructive techniques presented in GC are general, the results 
in GC are slanted toward single-level, univariate applications of covariance modeling. In 
accordance with terminology adopted in GC, and to qualify the extension of single-level 
univariate covariance functions to multiple levels, the nonseparable covariance functions 
constructed in this article will be referred to as multi-level. Multi-level covariance func- 
tions were introduced in GC, Secs. 3(c) and 4, where multi-level extensions of single-level 
covariance functions are described using graphs and analytic formulas valid for special 
cases. This article both extends and complements GC, with greater emphasis given below 
on relations between the three-dimensional Fourier spectrum and the covariance functions 
in physical space. A new, first order autoregressive covariance (FOAR) is developed as 
an example of a continuous, nondifierentiable multi-level covariance. The third-order au- 
toregressive covariance (TOAR) obtained by self-convolution in GC Sec. 4(b) is rederived 
using Fourier inversion, yielding formulas that more clearly show the relation between 
the multi-level TOAR and FOAR covariances. The multi-level piecewise rational covari- 
ance developed in GC Sec. 4(c) is discussed, with closed form expressions for the most 
general form of this covariance given in the Appendix. Schur products of single- and 
multi-level covariances are applied to yield additional multi-level covariances. Combining 
the multi-level covariance functions constructed below with single-level covariance func- 
tion examples constructed in GC yields many additional multi-level covariance functions. 
Multi-level powerlaw covariance functions are also constructed. The powerlaw covariances 
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have continuous derivatives of all orders, and thus are well-suited for applications requir- 
ing broad correlation functions. In the past several years, different forms of the multi-level 
powerlaw covariance have been applied to model both forecast and observation error co- 
variances within the PSAS. Section 3 describes the powerlaw covariance applied in the 
PSAS to produce tropically confined and longitudinally symmetric wind increments from 
spatially sparse wind observations in the tropical lower stratosphere. This new covariance 
model, which also has a large troposphere to stratosphere gradient in wind length scale, 
has been successfully applied to physically force a Quasi-Biennial Oscillation (QBO) from 
wind analysis increments in the Goddard Earth Observing System, Version 4 (GEOS- 
4), the operational data assimilation system developed at NASA/GSFC. This will be 
called the baseline GEOS-4 system in what follows. Section 3 below discusses the QBO, 
describes the new PSAS covariance model developed to force it from sparse wind obser- 
vations, and also shows a significantly improved QBO representation when compared to 
the QBO obtained using the baseline GEOS-4 system. 

2 Construction of multi-level covariance functions: theory and 

EXAMPLES 

First we introduce two definitions that provide equivalent tests to determine covari- 
ance functions on any set W. Definition 2.1 below extends GC Definition 2.1 to the case 
of vector-valued random fields. 

Definition 2.1 The multi-level covariance function of a random vector field 
X(r) := [X 1 (r),X 2 (r),...,Y m (r)] T , tEW, 
is the rn x m matrix function 


F(r, s) :=<Y(r)Y(s) T >, 


r, s 6 W, 


( 1 ) 
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where 

Y (r) := [ li(r), F 2 (r), . . . , Y m ( r) ] T , F*(r) := X k (r) - < X k (r) >, (2) 

where X k (r), k = 1, 2, . . . , m, are real-valued random fields on W, and where < ■ > de- 
notes mathematical expectation. □ 

For continuous covariance functions, Definition 2.1 is equivalent to the test described 
below in Definition 2.2 that is based on grid evaluations (cf. Loeve 1963, pp. 466-467). 

Definition 2.2 The m x m matrix of functions 

{B h (t,s)}, r.seW, k,l = l,2,...,m, (3) 

is a multi-level covariance function on W if for each positive integer n, and for each 
choice of points ri,r 2 , - - . , r n in W, the mn x mn matrix 

{Buiri, rj)}, *, j = 1, 2, . . . , n, k,l = l,2,...,m, (4) 

is positive semidefinite. □ 

The following three integration formulas are needed for the construction of multi- 
level covariance functions given below. 


Lemma 2.3 The Fourier transform (on three-dimensional Euclidean space R 3 ) off( r) := 
where fo(r) = exp(— |r|/L)/|r|, is given by 


/(w) := 


J exp (— i 


w T r) 


exp(-i|r||/X) 


dr = 


47 xL 2 

1 + L 2 ||* 


l|2 * 


( 5 ) 


Proof: Applying the Hankel transform (cf. GC Eq. (2.31)) yields 


OC 

a 47 r f 

/o(||w||) = J xfo(x) sin(l|w||x) dx = 


47 rL 2 


1 + L 2 ||w| 


| 2 ‘ 


( 6 ) 
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Lemma 2.4 The Fourier transform (on R 3 ) of B\{t) = exp(— ||rj|/L) is given by 


Bi{w) = 


(l + X2|| W || 2 ) 5 


Proof: Apply GC Eq. (2.32): 


where Bq(w) = 2L(l + w 2 L 2 ) 1 is given by GC Eq. (4.3). Alternatively, use Folland 
(1992, p. 247, problem 2). □ 


Lemma 2.5 For a > 0, a > 0, b > 0, 


OO 

/ ( M,a) := j = - JL_ (exp(-a6) — exp(— aa)) , a #6, 

— OO 

oo 

r/ , f x sin(ax) dx na , . 

(a,a ’ a) := J = 2a eXp( “ aa) ' (8) 

— OO 

Proof: Both integrals can be evaluated by direct application of the residue theorem (cf. 
Churchill 1984, Sec. 60). Alternatively, to find I (a, b, a) for a yt b, differentiate both sides 


of the identity (cf. Churchill 1984, p. 175, problem 4) 


cos(ax) 

( x 2 + a 2 ) (x 2 + h 2 ) 


7r fexp(-ab) exp(— aa)\ 

dx= J-Tv( i - a J 


with respect to a. The expression for I (a, a, a) can be found by taking the limit of both 
sides of /(a, b, a) as b tends to a. □ 

Example 2.6 below extends the so-called powerlaw of GC Sec. 4(d). These func- 
tions are everywhere infinitely differentiable and thus axe useful in applications where 
smoothness is required. In Section 3 below, a multi-level powerlaw function with long 
stratospheric length scales is shown to produce the approximately zonal analysis incre- 
ments for wind needed as forcing for a zonally symmetric quasi-biennial oscillation (QBO) 
from observational wind data. Multi-level covariance functions with zero, two, and three 
derivatives at zero separation are constructed following the multi-level powerlaw example. 
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Example 2.6 Given any length scales Lk > 0, k = 1, 2, . . . , m, define 

Pki (r, s) := ' - - a , t,s£R 3 , L kl := 0.5(L fc + Lj). 

Ljti 1 + 0.5(||r - s||/Ih) 

I7ie m x m matrix of Junctions 

{Pki( r,s)}, t,s£R k,l = 1,2, ...m, 
is a multi-level (powerlaw) covariance function. 


Proof: To verify Definition 2.2, we must show that 

m n 

£ £ P k l{Ti,Tj) C ik Cjl > 0 
fc,/=l i,j=l 

for arbitrary constants c,*, i = 1, 2, . . . , n, k=l,2, . . . ,m. Applying Lemma 2.3 
L = 1 (interchange the roles of r and w) yields 




1+IM 


| 2 - 


Using Eq. (13) and the change of variable (L k i\/2) u = w yields 


2TTPkl(Ti,Tj)= L k L t 


J exp(- 


i u T (ri - tj)) 


exp(-L w v / 2|]u||) 


du. 


lull 


Using the identity 

m n 

l cTz l ~ ^2 LkLt exp(-\\u\\L kl V2) expi-i^iri-Tj^CikCji, 

fc,i= 1 i,j=l 

where the mn-vectors c and z are defined by: 


C - [Cll ? ^21 j • • * j Cn 1 1 Cl2j * * ■ j ^n2 j • • • ) Ci jtj, • . . , Ctjitj.] , 

2 • [^lli -^21 1 - * - 5 ^filj ■2'12 i * * ■ 5 • * * ? • ■ * ? i 

Zji := L/ exp(m T rj) exp(-0.5 ||u||I,j\/2), j = 1, . . . , n, l = 1, . . . , m. 


in Eq. (14) yields 

m n , I T i 2 

^ ^ ^ 2| j 

2 " z2 z2 P u( T i, T j) CikCjl = / -jj-du > 0. 
k,l = 1 i,j = 1 ll U l‘ 


( 10 ) 


( 11 ) 


( 12 ) 

with 


(13) 


(14) 


(15) 


(16) 


( 17 ) 
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This completes the proof. □ 

Theorem 2.7 below generalizes GC Theorem 2.10 to a necessary and sufficient con- 
dition that characterizes all homogeneous multi-level covariance functions that are in- 
tegrable over R N . This result was first proven by Cramer (1940; cf. Yaglom 1987, 
pp. 314-315). A special case of this result can be applied to prove that the examples 
constructed below using cross-convolutions over R 3 are multi-level covariances- see the 
remarks following Thm. 2.7. 

Theorem 2.7 Define 

{B fc j(r,s)}, r ,s&R N , k, l = 1, 2, . . . , m, (18) 

where each Bm(t,s) is a homogeneous function on R N , and where each Bki (r, 0) is 
integrable over R N . Then {Bti{ r, s)} is a multi-level covariance function if, and only if, 
the m x m spectral density matrices 

J {Bkt( r, 0) exp(-t w T r) }dr, k, / = 1, 2, . . . , m, (19) 

are positive semidefinite for each w € R N . □ 

Remarks: Theorem 2.7 and GC Theorem 2.9 together imply that the m x m matrix of 
convolutions over R N : 

{(hk*hi)( r-s)}, r, s eR N , k, l = 1, 2, . . . , m, (20) 

is a multi-level covariance function, if each A* is a real- valued, radially symmetric function 
in L 1 (R N ) f! L 2 (R n ). This follows from the fact that the m x m spectral density matrix 
of Eq. (20) is the rank-one outer product 

{A fc (w)Aj(w)}, k,l = 1,2, ... ,m, (21) 

(cf. GC Theorem 2.9) and thus is positive semidefinite for each w € R N . 
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Observe that Thm. 2.7 was not applied in the proof of Example 2.6 above, since the 
multi-level power law is not integrablc over R 3 . □ 

The convolution integrals in Eq. (20) for the case of R 3 were reduced to a simpler 
form in GC Sec. 3(c). Several worked examples where these convolution integrals are then 
evaluated in closed form were given in GC Sec. 4. Three examples are provided below to 
complement and extend results presented in GC Sec. 4. The multi-level FOAR covariance 
given in Example 2.8 is new. It is non-differentiable at the origin, having cusp-like behav- 
ior, and thus is useful for modeling covariances with a slowly decaying spectral density 
(cf. Menard et al. 2000, Sec. 5). Compactly supported multi-level TOAR-like covariances 
were constructed in GC Sec. 4(b) through self-convolution of FOAR-like functions over 
R 3 , however, general formulas were not given in the case where the length scales are dif- 
ferent. Example 2.9 below provides explicit, simplified expressions for the multi-level 
TOAR covariance considered in GC Sec. 4(b). Recall that the TOAR covariance has 
four continuous derivatives at the origin. Multi-level compactly supported 5th-order 
piecewise rational functions extending GC Sec. 4(c) are given in Example 2.10 below. 
These functions have two continuous derivatives at the origin. 


Example 2.8 Define 


r* /_ ^ 2\/Z?V /expH|r - S \\/L k ) - exp(— ||r - s||/L,) \ _ , _ r /r 

F " (r ’ S) ' (U? - Lf) { jir^sji 


Fulr. r) := x5E, L„ := 0.5(1* + £,). 

Dkl 

Fkk (r, s) := exp(-||r-s||/Ljfe), k, l = 1 , 2, . . . , m. 


( 22 ) 


The m x m matrix of functions 


{-F*z(r,s)}, r, sGR 3 , k, l = 1, 2, . . . m, 


(23) 


is a multi-level (FOAR) covariance function. 
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Proof: To apply Thm. 2.7, it is shown that the multi-level FOAR function is a matrix 
of cross-convolutions over R 3 . To this end , define 




!|r||v2^It 


(24) 


and apply Lemma 2.3 to obtain 


A(w) = 


y/ 8irL k 3 

i+x* a iiwir 


(25) 


Lemma 2.4 shows that the Fourier transform of F kk (r, 0) is the square of f k ( w): 

8-rrL k 3 


/ F kk ( r, 0) exp(-iw T r)dr = — 

(‘ 


(26) 


+ V||w|| 2 ) 

so that F kk ( r, s) = (f k * f k )( r — s). To complete the proof, it is shown that F k i( r, s) is 
the inverse Fourier transform over R 3 of /*( w) jj(w), i.e., 

F k i (r, s) = (/* * /i)(r - s) = ^ J /*( w) / z (w) exp(iw T (r - s)) dw. (27) 

Let Ji be the value of the integral in Eq. (27), and let a k := 1/L k , at ~ 1/Li, where 
a k ^ai- Applying the Hankel transform yields 

OO 

2 [ x sin(||r — s||x) 


h = 


/ 


7rj|r — s\\y/LkLi J (x 2 + ai 2 ) (x 2 + a* 2 ) 

— OO 

/(a/,a fc ,||r-s[|) 


dx 


7r||r - s\\y/L k Lt 

2 ( exp(— ||r - s|la fc ) - exp(-||r - s||a z ) 


||r-s||vOTV 

= F ki ( r,s), 


a ( 2 - a* 2 


(28) 


where Lemma 2.5 was applied to evaluate the integral. □ 


Example 2.9 Define 


T k i(r, s) ~ -— 2 Lk2L [ 2 (. L k F kk (r , s) + L t F a { r,s) - 2 v^F M (r,s)), L* # L,, 

(it - L z ) v 


Tfcjfe (r, s) 1 + + 




3 V 


F kk (r, s), fc = l,2,...,m. 


(29) 
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The m x m matrix of functions 


(T w (r,s)}, r, s € f? 3 , k, l = 1, 2, . . . m, 


(30) 


is a multi-level (TOAR) covariance function. 


Proof: To apply Thm. 2.7, it is shown that the multi-level TOAR function is a matrix 
of cross-convolutions over R 3 . To this end, define 


«,(>) := * = 1,2,...,™, 




(31) 


and apply Lemma 2.4 to obtain 


5Jfe(w) = 


&V ttL/c 


(l +L fc 2 i|w|| 2 ) 


(32) 


First it is shown that 


Tfcfc(w) = 5t(w) 2 = 


64x Lk 


64x0* 5 


(33) 


(l+L* 2 ||w|| 2 ) (o* 2 + !l w || 2 ) 

where a k := 1/L*, implying that X**(r, s) = (<?* * <?*)( r - s). Applying the Hankel trans- 
form yields 

( 9k*gk)(r-s ) = ^ J 5*(w) 2 exp(iw T (r-s))dw 

o° 

16a* 5 f x sin(||r — s||x) 

x||r-s|| J ( x 2 

OO 

/ 


8a* 5 f cos(||r — s||x) 


{x 2 + a* 2 ) 4 
3 


dx 


(34) 


3 x y (*2 + ajt 2) = 

— oo 

where the third equality in Eq. (34) was obtained using integration by parts. To simplify 
Eq. (34), recall the relation between the one-dimensional SOAR function and its spectral 
density function (cf. Daley 1991; Churchill 1984, p. 176, problem 6): 

OO 

h{a k , ||r||) •- J = ^3 (* + a *H r ll) exp(- Q *||r||). (35) 
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Differentiating both sides of Eq. (35) with respect to a* and substituting into Eq. (34) 
yields 

*2.(1 ^ 

(9k*9k)( r-s) = — ^-h'{a k , i|r-s||) = T kk { r,s). (36) 

Now let a-i := 1/E/, with a k ^ a; . To complete the proof, it is shown that 

Tu( w) = P/fe(w)p/(w) = ^^ 2 ? L ' 71 

(l + Efc 2 |Jw|| 2 ) (l + i, 2 j|w|| 2 ) 

647r>/a fc 5 ai 5 

/ 9\ 2 / 9\ 2 ' 

( a * 2 + IMI ) (a / 2 + ||w|| 2 J 

Applying the Hankel transform yields 


647ry/Et 3 E/ 3 

(l + E fc 2 ||w|| 2 ) 2 (l + E, 2 i|w|| 2 ) 

647r>/afc 5 ai 5 


(9k*9i)( r-s) = / </k(w)$j(w) exp(iw T (r - s)) dw 

OO 

_ 16yafc 5 ai 5 r xsin(||r — s||x) ^ 

7 r ll r — s ll J (x 2 + a* 2 ) 2 (x 2 + a / 2 ) 2 

— OG 

= li — ^2 (f(ofe,o fc ) +J(oi,a/) - 2/(a fc ,a/)), (38) 
fliF — s||(o* 2 — a/ 2 ) v > 

where the third equality is obtained by expanding (x 2 + a fc 2 )~ 1 (x 2 + a/ 2 ) -1 into partial 
fractions, squaring the result, and applying Lemma 2.5, where I (a, b) := /(a, 6, ||r - s||) 


is abbreviated. Substituting 


V L k Li F k i (r, s) = 


7r r-s 


I(a k ,at, ||r-s||) 


(39) 


into Eq. (38) and simplifying yields the result. □ 


Example 2.10 Given scalars a k and c, define 

( (2 (a* - l)!|r||/c+ l)n fe 0 < !|r|| < c/2 
hfc(r) := < 2a*(l - ||r||/c)n* c/2<|jr||<c, k = 1, 2, . . . , m, (40) 

l 0 c < ||r|( 

where n k := (2 + 44a| + 6a*) 1/,S . Closed form expressions for the convolution integrals 


h k * hi over i? 3 


{h k *hi){ r-s) = / 


' /i(l|r-s||/c)n*n/, 
/ 2 (||r — s\\/c)n k ni, 
MW T ~ s \\/c)n k ni, 
/4(l|r- s|| /c)n k ni, 
0 


0 < )|r — s|| < c/2 
c/2 < ||r — s|| < c 
c < ||r — s|| < 3c/2 
3c/2 < ||r - s|| < 2c 
2c < llr -sll 


(41) 
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are given in the Appendix. The case a* = a/ is illustrated in GC Figs. 7,8. If a* = a/ = 
1/2 then Ea. (41) reduces to GC Eq. (4,10). □ 

Theorem 2.11 below provides an important practical technique for the construction 
of multi-level covariance functions using the Schur product. 

Theorem 2.11 IfC{ r, s ) is a covariance function on a set W and 

{B kl { r,s)}, r, s€ W, k,l = l,2,. . . ,m, (42) 

is a multi-level covariance function on W, then 

{C(r,s) Bju(r,s)}, r.sef, k, l = 1, 2, . . . , m, (43) 

is a multi-level covariance function on W. 

Proof: Since C(r, s) is a covariance function on W, the mn x mn matrix 

{Cfcl(r,s)}, Cu :=C, r, s € W, h, l = 1 , 2, . . . , m, (44) 

is a multi-level covariance function on W: 

m n n m 

Y Y C(T l ,r J )b lk b J i = Y c ( T i’ r i) Y hikb i l ( 45 ) 

k,l= 1 i, j—1 i,j= 1 fc,Z=l 

n tjx m 

— yi bji 

i,j= 1 k— 1 l~ 1 

n 77i 

= Y C(ri, rj) fcPj > 0, Pi~Y bik ' 

t,j=l k=l 

The theorem follows by applying the Schur product theorem to the two positive semidef- 
inite mn x mn matrices {Cfc;(r ; , rj)} and {B k i (r ; , rj)}. □ 

Example 2.6 and Theorem 2.11 above provide tools needed to construct nonseparable 
extensions using Schur products. One such extension, summarized in Example 2.12 below, 
has been applied in the Physical-space Statistical Analysis System (PSAS) (Cohn et al. 
1998) for the past several years. 
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Example 2.12 The m x m matrix of functions 

P fc /(r,s) C(r, s), r, s £ R 3 , k, l = 1, 2, . . . , m, (46) 

where Pki (r, s) is the powerlaw defined in Example 2.6 above, and where C(r, s) is any 
correlation function on R 3 , is a multi-level covariance function. 

Remark: The PSAS application of Example 2.12 in the baseline GEOS-4 system uses 
C represented by GC, Eq. (4.10). □ 

Theorem 2.13 below generalizes GC Theorem 2.11 to a necessary and sufficient con- 
dition that characterizes all continuous isotropic multi-level covariance functions on the 
(N — l)-dimensional sphere S N ~ 1 . This result was first proven by Hannan (1970; cf. Ya- 
glom p. 387). 

Theorem 2.13 Consider the m x m matrix of functions 

(J3 w (r,s)}, r.seS* -1 , k, l = 1, 2, . . . , m, (47) 

where each Bki{r , s) is a continuous isotropic function on S N ~ l . Then { B k i(j , s)} is a 
multi-level covariance function if, and only if, 

OO 

{Bki (r, s)} = Y, A iCi^- 2 ^ 2 (r T s), (48) 

t=0 

where each A* is a positive semidefinite m x m matrix, and where C^ N ~ 2 ^ 2 are Gegen- 
bauer polynomials of degree m and order (N — 2)/2. □ 

Remark: Isotropic multi-level covariance functions on S N ~ 1 are readily obtained using 
the above results developed for R N simply by restricting the points from R N to S^" 1 . 
This point is further elaborated in GC Sec. 2(c). □ 
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3 AN APPLICATION: TROPICAL STRATOSPHERIC WINDS 

The QBO dominates the dynamics of the equatorial stratosphere between about 
70 hPa and 10 hPa, with downward propagating, zonally symmetric easterly and westerly 
wind regimes. Observations of the QBO indicate that it is symmetric about the Equator 
with a maximum amplitude of about 20 ms -1 and has an approximate Gaussian half- 
width of about 12° in latitude (cf. Holton 1992, pp. 424-425). Dynamical aspects of 
the QBO are generally illustrated using a time-height cross section of monthly-mean 
equatorial zonal wind (e.g. Baldwin et al. 2001). It is commonly accepted that the QBO 
is approximately zonally symmetric (Belmont and Dartt 1968), although a recent study 
by Hamilton et al. (2004) showed that zonal asymmetries can be found at 10 hPa during 
the westerly phase of the QBO. This point will be further elaborated below. 

Efforts to understand the QBO using observational data are complicated by the lack 
of available in-situ wind observations in the near-equatorial stratosphere above 50 hPa. 
Radiosonde reports are scarce, since many balloons burst near the very cold tropical 
tropopause (Hamilton et al. 2004). As a result, there have been difficulties in generating 
accurate QBO signals in some assimilated datasets (cf. Randel et al. 2004). A common 
problem is that the amplitude of the QBO is often underestimated in comparison with 
radiosonde measurements (Pawsonand Fiorino, 1998) although the ECMWF Reanalysis- 
40 (ERA-40) does very well (cf. Randel et al. 2004). 

Figure 1 illustrates that radiosonde observations available to reproduce a QBO are 
concentrated in the Fair East. The top panel of Fig. 1 shows the station averaged value of 
the 10 hPa zonal wind, at 00 GMT and rounded to the nearest integer, for each station 
reporting 5 or more times in November 1991. The influence these observations have on the 
assimilated zonal wind is discussed below. The middle panel shows the number of reports 
corresponding to each location shown in the top panel. For instance, the average wind 
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speed of 7 ms 1 in the top panel at the Singapore radiosonde (1.37° N, 103.98° E) is over 
26 of 30 possible 00 GMT reports. The bottom panel shows the number of radiosonde 
stations reporting 5 or more times at 50 hPa. A typical distribution of radiosonde reports 
on mandatory pressure levels between 10 hPa and 50 hPa is less dense than the lower 
panel, and more dense than in the middle panel. Dynamical properties of the QBO are 
often illustrated using time series of the zonal wind observed at the Singapore radiosonde, 
due to the sparsity of available observations, and since the reports from this station are 
both consistent and of high quality. A centered 10-day moving average of zonal wind 
observations from 10 hPa to 50 hPa at Singapore is shown in Fig. 2 from May 1991 
through November 1995. Several observed properties of the QBO are clear in this figure. 
For example, easterly and westerly wind regimes alternate, the maximum amplitude is 
approximately 20 ms -1 from 10 hPa to 30 hPa, and the downward propagation of zonal 
wind occurs without loss of amplitude from 10 hPa to 30 hPa. 

To generate a QBO in GEOS-4, tropically confined and longitudinally elongated 
zonal wind analysis increments are needed to effect the rapid QBO phase transitions 
illustrated in Fig. 2. Two multi-level covariance models were tested in GEOS-4. The 
model used to improve the QBO representation, henceforth known as model A, is a 
powerlaw in the form of Example 2.6 with tropically confined and longitudinally elon- 
gated wind correlation length scales. Model B is the so-called windowed powerlaw used 
in the baseline GEOS-4 data assimilation system, in the form given in Example 2.12. 
A static PSAS analysis using a single Singapore radiosonde sounding at 1200 GMT on 
21 September 1999 is described below for the purpose of illustrating key differences be- 
tween models A and B. Figure 3 shows the correlation length scale differences between 
models A and B (left panel) and the observed-minus-forecast (OMF) residuals used in 
this static PSAS analysis (right panel). The length scales are identical for p > 200 hPa. 
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Model A has much larger length scales than model B for p < 200 hPa. The PSAS zonal 
wind analysis increments produced using models A and B are shown in Figs, 4 and 5. At 
2° N (Fig. 4), the increments are similar where the length scale profiles shown in Fig. 3 
are identical, but their longitudinal structure reflects the stratospheric length scale differ- 
ences. The analysis increments near Singapore are, as expected, a fraction of the nearby 
OMF values (Fig. 3). Such increments illustrate a typical action of the gain matrix on 
the OMF vector (cf. Daley 1991, p. 377). Figure 5 shows horizontal maps at 10 hPa of 
the increments seen in Fig. 4, illustrating the approximate longitudinal symmetry of the 
analysis increments produced using model A, and the more confined meridional scale. 

Assimilation experiments A and B, using GEOS-4 with covariance models A and 
B, respectively, were configured to study the impacts on the assimilated QBO of using 
the elongated correlation length fields in the tropical stratosphere. The analysis focuses 
on the representation of the winds and shear zones at the onset of the westerly wind 
regime at 10 hPa, between November 1991 and January 1992. The two experiments were 
initialized with identical states on 1 May 1991. To ensure that both experiments used 
the same tropical wind observations, the GEOS-4 quality control settings were updated 
to accept all radiosondes in the tropical lower stratosphere. 

The first major differences between the two experiments occur at the onset of the 
10 hPa phase transition in late October 1991 (see Fig. 2). Figure 6 shows aspects of the 
zonal wind at 00 GMT on 30 November 1991, approximately one month following the 
onset of westerly winds at Singapore. The upper two panels illustrate typical differences 
between the vertical structure of the zonal- mean zonal winds in the two assimilation 
experiments. Experiment A represents a plausible zonal-mean wind field with greater 
sensitivity to the observed vertical shear (Fig. 2) than Experiment B. The strong easterlies 
in Experiment B resemble typical observations seen several months prior to the QBO 


18 


Gregory Gaspari, Stephen E. Cohn, Jing Guo, and Steven Pawson 


phase transition. The difference between the two zonal averages reveals the substantial 
positive impact obtained by stretching the wind correlation length field in longitude. 
The shear zones in Experiment A seen near 10 hPa and 60 hPa reflect the wind shear 
present in the observations on 30 November 1991. The large differences near 10 hPa 
illustrate a systematic difference between the two experiments during the entire QBO 
phase change, while the differences near 60 hPa are primarily due to less systematic 
variations in the radiosonde observations. The profiles shown in the lower right panel of 
Fig. 6 illustrate differences between the local analysis near Singapore and the zonally 
averaged zonal-mean wind field in both experiments. In Experiment A, there is close 
agreement between the local analysis near Singapore and the zonal average on the latitude 
circle. The corresponding profiles for Experiment B are less s imil ar due to the short 
correlation length scales. Lack of sufficient “spreading” of the analysis increments in 
Experiment B causes the wind shear near 10 hPa to be misrepresented. Figure 7 shows the 
10 hPa winds and their analysis increments for Experiments A and B, averaged over the 
00 GMT data for November 1991. This figure clearly shows the introduction of westerly 
winds near the Equator in Experiment A, and the remains of a strong easterly flow in 
Experiment B. There is some zonal asymmetry in Experiment A, with two wind maxima 
of 12 ms -1 over the Southern Indian Ocean and 9 ms -1 over the eastern Pacific, with 
weak easterlies over the Atlantic. The analysis increments shown in the lower panels can 
be viewed as typical corrections applied every 12 hours (when radiosondes are present) to 
the background zonal wind. The analysis increments in Experiment A are small and quite 
symmetric near the equator. The small size of the increments (fractions of one meter 
per second) indicates the model and observations are working in harmony. Imposing 
the zonally elongated analysis increments means that initial changes are spread over 
a wider longitude range, allowing the model to sustain a slowly increasing region of 
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westerlies. This is not the case in Experiment B, where the localized analysis increments 
apparently have no lasting effect on the zonally averaged zonal wind field. Since these 
analysis increments are both localized and intermittent (every 12 hours), advection of 
the strong easterlies quickly leads to their local replenishment, so that a large difference 
between model forecast and observations remains and the model is not able to develop 
the observed westerly winds. By January 1992, the transition to westerlies in Experiment 
A has occurred at all longitudes (Fig. 8), with a displacement of the westerly maxima to 
the North of the Equator. In experiment B, the wind remains easterly and the analysis 
increments retain a complex structure. 

4 Concluding Remarks 

Theoretical and practical mathematical results were developed to construct covari- 
ance functions dependent on variable length fields. These general results could potentially 
be applied in many application areas requiring parametrized covariance models having 
a modest degree of additional flexibility beyond homogeneous and isotropic models. The 
specific application to improve the QBO illustrates the ability of the covariance model 
to incorporate both tropospheric and stratospheric length scales appropriate to different 
regions of the atmosphere. Other applications will be explored in future work. 

5 Appendix 

The functions /i, . . - , fi defined in Eq. (41) are given by: 

/i(z) = bisz 5 + &i4Z 4 + bi^z 3 + 6i2Z 2 + bio, (49) 

/z( z ) = b25Z: 5 + 624-2 4 + &23~ 3 + b22^ 2 + &2lZ + 620 + d2lZ 1 , 
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f3( z ) — &35Z 5 + &34Z 4 + bzsZ 3 + & 32 Z 2 + 631 Z + 630 + d^iZ X , 

fi(z) = (2z 5 - 12z 4 + 15 z 3 + 40 z 2 - 120z + 96 - IGz -1 ) , 

bis 16(— 3 — 7a*a/ + 4a* + 4a f )/3, 

614 := 16(1 + 2 a*a/ — a* — a/), 

&13 := 10(1 + 8 a*a/ — 2 a* — 2 a/), 

612 := -40(1 + 8 a*a/ - a* - a/)/3, 

610 *— 2 -t- 44a* a/ 4- 3a* -P 3a/, 

625 := 16(1 + 5a*a/ - 3a* - 3a/)/3, 

624 ;= 8 (— 2 — 8 a*a/ + 5a* + 5a/), 

623 := 10, 

622 := 20(2 + 20a*a/ — 11a* — lla/)/3, 

&21 := 5(— 4 - 36a*a/ + 13a* 4- 13a/), 

620 := (16 + 204a*a/ - 35a* — 35a/)/2, 

£^21 := (—8 — 84 a*a/ + 29a* + 29a/)/12, 

635 := 16(— 3a*a/ + a* + a/)/3, 

634 := 8 ( 8 a* a/ - 3a* - 3a/), 

& 33 := — 20 ( 2 a*a ; - a* - a/), 

632 — 20(20a*a/ — 9a* — 9a/)/3, 

631 := 5(44a*a/ - 27 a* - 27a/), 

630 := — (244a*a/ — 189a* — 189a/)/2, 

<^31 : = 


(460a*a/ — 243a* — 243a/)/12. 


(50) 
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Figure 1. The top panel is average value of the 10 hPa zonal wind over Nov. 1991 at 0000 GMT, for all 
tropical radiosonde stations reporting 5 or more times in Nov. 1991. The 9 ms -1 value at 3.75° N, 102° E 
is the reporting frequency weighted average over the two stations at (3.1°, 101.65°) and (3.62°, 103.22°). 
The location (3.75°, 102°) is the center of the 2.5° x 4° lat-lon box used for the averaging. All other 
values are shown at the reporting station locations. The second panel is the total number of 10 hPa 
zonal wind reports at 00 GMT over Nov. 1991. The combined number of reports from (3.1°, 101.65°) 
and (3.62°, 103.22°) is 12, shown at (3.75°, 102°). The third (lowest) panel marks the location of all 
tropical radiosonde stations with 5 or more 50 hPa wind reports at 0000 GMT in Nov. 1991. Note the 
negative averages of -13 ms -1 at (7.33°, 134.48°) and -10 ms -1 at (—14.33°, 189.28°). 
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Figure 2. Ten day centered moving average of zonal wind radiosonde observations in ms -1 at the 

Singapore radiosonde. 
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Figure 3. The left panel is the streamfunction length scales in models A and B, where the longer length 
scales are used in model A. The x-axis is the length scale in kilometers and the y-axis is pressure in hPa. 
The right panel is the observed-minus-forecast zonal wind in ms -1 (x-axis) at Singapore used in a static 
PSAS analysis to produce the analysis increments shown in Figs. 4 and 5. 
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Figure 4. Zonal wind analysis increments at latitude 2° N in ms -1 from a static PSAS analysis using 
a single radiosonde profile at Singapore with observed-minus-forecast zonal winds shown in Fig. 3. The 
upper panel shows the increments due to model A. The lower panel shows the increments due to model 

B. The contour interval is .5 ms -1 . 







panel shows the increments due to model A. The lower panel shows the increments due to model B. The 


contour interval is .2 ms 1 . 
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(d) Profiles 


latitude 


Figure 6. The upper panels show the longitudinally averaged zonal wind field in ms - 1 at 00 GMT on 
30 November 1991. Experiment A is shown on the left and Experiment B is shown on the right. The 
y-axis is pressure in hPa. the x-axis is latitude, and the contour interval is 10 ms -1 . The difference of the 
upper left and right panels is shown in the lower left panel with a contour interval of 5 ms -1 . Profiles 
of the zonal wind field in ms -1 near Singapore are shown in the low'er right panel, where the x-axis is 
the zonal wind and the y-axis is pressure in hPa. The dotted (dashed) solid curve is the zonal wind for 
Exp. A (Exp. B) at latitude 1° X and longitude 103.75“ E. The thick (thin) curve is the longitudinally 
averaged zonal wind for Exp. A (Exp. B) at latitude 1” X. 
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Figure 7. The monthly averaged zonal wind for Exp. A and B in ms -1 at 10 hPa at 00 GMT over 
November 1991 is shown in the upper two panels. The monthly averaged zonal wind analysis increments 
in ms -1 for Exp. A and B at 10 hPa at 00 GMT over November 1991 is shown in the lower two panels. 
The averages for Exp. A (Exp. B) are shown in the left (right) panels. The contour interval for the upper 

(lower) panels is 3 ms -1 (.2 ms -1 ). 
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Figure 8. As in Fig. 7, for January 1992. 




Popular Summary 

Construction of Covariance Functions with Variable Length Fields 

Gregory Gaspari, Stephen E. Cohn, Jing Guo, Steven Pawson 

Geophysical data assimilation involves using principles of estimation theory to merge diverse 
observations and short-term model predictions to yield so-called "optimar analyses of the system 
state. The optimization in assimilation is effected using the best possible weighting of model errors 
against observation errors. The ability to accurately determine this weighting is arguably one of the 
least well understood aspects of the assimilation, due primarily to the diverse range of conditions 
encountered in the atmosphere. For this reason, the quality of an assimilation is typically evaluated 
experimentally, based on the geophysical response to perturbations to one or more of the following: 
the statistical analysis, the observing system, and the model. The former two components are 
interrelated in the sense that covariance modeling assumptions determine the spatio-temporal 
influence of observations through so-called "analysis increments" that act to iteratively correct the 
model state over the course of the assimilation. In meteorological applications, analysis increments 
are typically produced using background and observation error covariance models applied to 
determine the spatio-temporal influence of a diverse set of in-situ, or "point" measurements (e.g., 
balloon soundings, aircraft observations) with a diverse set of satellite observations (e.g., space- 
based radiance data). Effective analysis increments should cause high-quality observations to be 
distributed over length scales appropriate to the region of the atmosphere, and appropriate for the 
spatio-temporal data coverage. Geophysical length scales differ widely between the atmospheric 
boundary layer, the free troposphere and the stratosphere. A practical determination of the 
correlation length scale should account for these differences, while also accounting for differences 
between data-dense and data-poor regions of the atmosphere. 

This paper presents new theoretical developments for a formal mechanism of incorporating 
geographically- and height-dependent length scales in the assimilation, extending earlier theoretical 
developments of Gaspari and Cohn (1999). 

An application is presented for the assimilation of tropical wind observations, where the length 
scale in the sparsely observed stratosphere is much longer th an that in the more densely observed 
troposphere. It is demonstrated, using the GMAO's GEOS-4 data assimilation system, that using a 
zonally extended length scale in the tropical stratosphere leads to a much more realistic quasi- 
biennial oscillation than with the shorter length scales used in the standard system. 



